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Ice sheets appeared in the northern hemisphere around 3 Ma (million years) ago and 
glacial-interglacial cycles have paced Earth's climate since then. Superimposed on these 
long glacial cycles comes an intricate pattern of millennial and sub-millennial variability, 
including Dansgaard-Oeschger and Heinrich events. There are numerous theories about 
these oscillations. Here, we review a number of them in order to draw a parallel between 
climatic concepts and dynamical system concepts, including, in particular, the relaxation 
oscillator, excitability, slow-fast dynamics and homoclinic orbits. Namely, almost all 
theories of ice ages reviewed here feature a phenomenon of synchronization between 
internal climate dynamics and astronomical forcing. However, these theories differ in 
their bifurcation structure and this has an effect on the way the ice age phenomenon 
could grow 3 Ma ago. All theories on rapid events reviewed here rely on the concept 
of a limit cycle excited by changes in the surface freshwater balance of the ocean. The 
article also reviews basic effects of stochastic fluctuations on these models, including the 
phenomenon of phase dispersion, shortening of the limit cycle and stochastic resonance. 
It concludes with a more personal statement about the potential for inference with simple 
stochastic dynamical systems in palaeoclimate science. 

Keywords: palaeoclimates; dynamical systems; limit cycle; ice ages; 
Dansgaard-Oeschger events 



1. Introduction 

The Pliocene and the Pleistocene cover approximately the past 5Myr. The 
climatic fluctuations that characterized this period may be reconstructed from 
numerous natural archives, including marine, continental and ice core records. 
These archives show a complex climate history. Ice sheets appeared in the 
northern hemisphere around 3-3.5 Ma (million years) ago [1,2]. The volume 
of these ice sheets fluctuated with the variations of the seasonal and spatial 
distributions of incoming solar radiation (insolation), which are induced by 
changes in the geometry of the Earth's orbit and the angle (obliquity) between 
*michel. cr uciflx(Q)uclouvain . be 
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Figure 1. Climatic fluctuations over the Late Pleistocene. Ice ages are reconstructed using the 
oxygen isotopic ratio of the calcite shells of benthic foraminifera [9]. Within the last ice ages, large 
temperature fluctuations were recorded in Greenland, here depicted by changes in the oxygen 
isotopic ratio of ice [10]. 



Earth's equator and the ecliptic [3-5]. This is called the astronomical forcing. 1 
Glacial cycles had an average duration of about 40000 years [6] until about 
800 000 years ago. The dominant period of glacial cycles increased around 800 000 
years ago and this is referred to as the Middle Pleistocene Transition. Data and 
models about the Middle Pleistocene Transition are reviewed by Clark et al. [7]. 
Time-series analyses based on band-pass filtering provide further evidence of the 
nonlinear nature of the climate response to the astronomical forcing, from about 
1.4 Ma ago [8]. The latest four glacial cycles, in particular, are distinguished by a 
pronounced saw-tooth time structure: ice accumulates over the continents during 
about 80000 years and then melts in about 10000 years (figure 1). 

Superimposed on these long glacial cycles comes a complex pattern of 
millennial and sub-millennial variability [11]. For example, the Greenland record 
features at least 20 events of abrupt rise and slower decline in oxygen isotopic ratio 
(a proxy for temperature) [12,13] and methane [14] during the latest glacial epoch. 
These events are known as Dansgaard-Oeschger events. They were found to occur 
from at least the last glacial inception [15] and the Antarctic ice core records 

lr rhe astronomical forcing will generally be taken into account here in the form of a normalized 
measure of insolation during the month or on the day of summer solstice at a northerly latitude, 
typically 60 or 65° N. This is a fairly complex, aperiodic signal, with dominant harmonics 
corresponding to the phenomena of precession (23 716, 22 428 and 18 976 years), and obliquity 
(41000 years) [5]. 
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provide evidence that they are characteristic of Pleistocene glacial climates [16]. 
Some of these events follow pulses of iceberg discharges into the North Atlantic 
Ocean, called Heinrich events [17-19]. Heinrich events and Dansgaard-Oeschger 
events have left climatic footprints all over the globe [20], including in Antarctica 
[16]. The current interglacial period is referred to as the Holocene. It is also 
characterized by millennial and centennial variability, mainly observed in the 
North Atlantic [21-23] , but of a much weaker amplitude than during the preceding 
glacial period. 

This paper reviews attempts to explain these fluctuations with concepts that 
originate in dynamical system theory. These are the concepts of limit cycle, 
synchronization and excitability. The central message of the paper is that current 
theories of ice ages and rapid events may often be interpreted in terms of generic 
deterministic models, which are also used in other areas of science such as biology 
and ecology. However, stochastic parametrizations are an essential part of any 
complex system model, and their effects on climatic oscillations have to be taken 
into account. 

Dynamical system theory entered palaeoclimate science with idealized models 
representing the response of ice sheets to the astronomical forcing. These models 
were directly derived from the physics of the ice sheet-atmosphere system [24-27] . 
Ghil & Childress [28], in particular, insisted on the interest of analysing such 
models in terms of bifurcation theory. For modelling the complex carbon cycle 
response, authors sometimes adopted a more heuristic approach by considering 
simple models and confronting the results of palaeoclimate evidence [29]. 

Nowadays, climate research is largely oriented towards large climate simulators 
(typically general circulation models), which are developed to include as many 
climate processes as possible. However, thinking in terms of dynamical system 
theory remains insightful. Indeed, the behaviour of a complex system at a 
certain spatio-temporal scale is in practice often dominated by a few leading 
modes, of which the dynamics may be captured fairly convincingly with a low- 
order dynamical system. Climate scientists are increasingly using this property. 
For example, they formulate simpler models to explain the seemingly complex 
behaviours observed in ocean-atmosphere simulators. Examples have been 
provided in recent years by focusing on interannual [30], centennial [31] and 
millennial [32,33] variability. In parallel, so-called hysteresis experiments, which 
aim at identifying the number of stable states in individual components of the 
climate system such as the ocean circulation [34] or ice sheets [35] , contribute to 
a dynamical system-founded understanding of the climate system. This approach 
may also help us to predict and communicate about the proximity of bifurcations, 
which may result in catastrophic climatic changes. Timmermann & Jin [36] 
termed our ability to anticipate bifurcation phenomena as predictability of the 
third kind, by reference to the predictabilities of the first and second kind 
originally introduced by Lorenz [37]. 

The article is structured as follows. Section 2 reviews some of the basic 
concepts of oscillator theory. This is no substitute for proper textbook reading, 
but the reader will find essential notions and definitions needed to understand the 
remainder. Section 3 reviews how these concepts enter theories of ice ages and 
rapid events. Section 4 discusses effects of stochastic fluctuations and, finally, 
§5 is a more personal statement about the potential for inference with simple 
stochastic dynamical systems in palaeoclimate science. 
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2. Vocabulary and elementary notions 

The reader will find an accessible introduction to dynamical system theory and 
concepts in Strogatz [38]. More formal background on oscillator theory, albeit a 
bit dated, is available in Guckenheimer & Holmes [39]. Bifurcation and oscillator 
theory is explicitly connected to climate theory in Ghil & Childress [28, in 
particular, ch. 12] and Saltzman [40, ch. 7]. Background on synchronization 
and an introduction to the phenomenon of excitability is available in Pikovski et 
al. [41]. Finally, the Scholarpaedia peer- reviewed website is an increasingly rich 
and authoritative source of information on dynamical systems. Only the notions 
essential for the present article are summarized here. 

Oscillator. The oscillator is a dynamical system that has a globally attracting 
limit cycle. In more simple terms, it oscillates even in the absence of an external 
drive. Here, we are interested in oscillators to describe climate phenomena, which 
involve dissipation of energy. The minimal model for a dissipative oscillator 
includes two ordinary differential equations, of which at least one is nonlinear. 

Relaxation oscillator. The relaxation oscillator is a particular kind of oscillator 
featuring an interplay between relaxation dynamics (generally fast) and a 
destabilization process (generally slow). The relaxation is the process by which 
the system is attracted to a region of the phase space. This evokes the relaxation 
of a spring. In a relaxation oscillator, the system continues to evolve slowly after 
the relaxation phase. During this slow evolution phase, the system's stability 
diminishes gradually until the system is ejected out of its relaxation state, either 
towards another relaxation state or to the same relaxation state via a dissipative 
loop. In this review, we will encounter three kinds of relaxation oscillators 
(figure 2): relaxation founded on slow-fast dynamics (involving a slow manifold); 
relaxations structured by a homoclinic orbit (involving only one relaxation state); 
and relaxations structured around a focus. More details are provided in the 
caption of figure 2. 

Excitability. An excitable system has a globally attracting fixed point (it does 
not oscillate spontaneously). However, an external perturbation may have the 
effect of exciting it. During this excitation, the system is being ejected far from 
its fixed point and then returns to it. 

Link between relaxation dynamics and excitability. In practice, it is often found 
that a relaxation oscillator may be transformed into an excitable system by a mere 
change in parameter, and vice versa. The reason is the following. A relaxation 
oscillation is often structured globally in the phase space; for example, by a slow 
manifold (figure 2a) or by one or several saddle points (figure 2b). Suppose now 
that the oscillation displayed by such a system ceases because a parameter has 
been changed. The system is then no longer an oscillator, but the 'backbone' of 
the oscillation dynamics are still latent in the phase space because the elements 
that structured the limit cycle (the slow manifold or the saddle points) have not 
disappeared. Consequently, the system may be run on a trajectory close to the 
defunct limit cycle if it is being pushed by some external force (the excitation) 
into the region of the phase space previously occupied by this limit cycle. This 
point is illustrated on the basis of slow-fast dynamics in figure 2d, but similar 
excitation dynamics generally occur near any kind of 'explosive bifurcation', that 
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Figure 2. Schematic of several forms of relaxation oscillations, (a) The vector space is structured 
by a slow manifold with several stable branches. All points of the state space are attracted towards 
the stable branches of the slow manifold (solid lines) along the fast direction, which is here the 
horizontal. The slow evolution consists of an upward or downward course along the slow manifold 
depending on whether the system lives on the right- or left-hand side of the null-cline of the slow 
variable. The relaxation oscillation consists of alternate jumps between the two branches of the 
slow manifold, (b) Trajectories are rapidly attracted towards a region of the phase space influenced 
by a saddle point. In the scenario displayed here, there exists a combination of parameters for 
which the limit cycle crosses the saddle point. In that case, the period of the orbit, which includes 
the saddle-node, is infinitely long. It is a 'homoclinic orbit', hence this particular bifurcation is 
named a homoclinic bifurcation. The homoclinic orbit only exists at the bifurcation point, but it 
influences the orbit when the parameter is close to the bifurcation. This is the reason why one refers 
to the 'ghost' of the homoclinic orbit. There is another scenario for which different saddle-nodes are 
connected to each other. The orbit and the associated bifurcation are then said to be 'heteroclinic'. 
(c) The relaxation oscillation is organized around a fixed point, with complex eigenvalues with a 
positive real part. The bifurcation giving rise to this orbit is a Hopf bifurcation, (d) One example 
of excitability, here depicted for a slow-fast system. The system resides in a stable space, but a 
fluctuation may cause an ejection out of the unstable (dashed lines) branch of the slow manifold. 
The system then loops all the way through the slow manifold before coming back to rest. (Online 
version in colour.) 
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is, bifurcations that give rise rapidly to a fully developed limit cycle. This includes 
homoclinic, heteroclinic and certain Hopf bifurcations (two examples follow and 
are illustrated in figure 6). 



3. Oscillators, relaxation and excitability in palaeoclimates 



(a) Models of ice ages 

(i) The Saltzman et al. models 

Saltzman established a theory in which ice ages are interpreted as a limit cycle 
synchronized on the astronomical forcing. Saltzman et al. [42] wrote a series of 
articles on the subject, starting with the introduction of the limit cycle idea 
[43] and synchronization hypothesis, the interpretation of the Middle Pleistocene 
Transition as a bifurcation [44], and the more complete models in the mid-1990s 
[45]. The full theory is developed in a book [40]. Here, we concentrate on two 
intermediate models [46,47]. They are called SM90 and SM91, by reference to the 
authors (Saltzman & Maasch) and the year of publication. The variables /, ill and 
6 are the continental ice mass, CO2 concentration and deep-ocean temperature, 
respectively. The reader is referred to the original publications for the meaning 
and value of the different parameters. They are not crucial here; it suffices to 
know that they are all positive. 
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In both models, the first equation describes the ice mass response to changes 
in CO2 (m) an d the astronomical forcing (F/(t)); Saltzman adopts the so-called 
Milankovitch view 2 that an increase in insolation causes a decrease in ice mass. 
Increases in CO2 or in ocean temperature have the same effects. 

The other two equations describe the dynamics of CO2 and the response 
of deep-ocean temperature to changes in ice volume. It is further assumed 
that the mean state of the climate varied slowly throughout the Pliocene- 
Pleistocene; in particular, in response to a 'tectonically driven' decline in the 

2 In fact, the view is introduced by Murphy [48] but it is developed mathematically in the 'canon 
of insolation' authored by Milankovitch [4], 
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Figure 3. Bifurcation diagrams of the Saltzman and Maasch models (a) SM90 and (b) SM91 as 
a function of F^, here treated as a constant control parameter. Note the difference in scales on 
both axes. Black lines are fixed points. Continental ice mass (/) is shown as a function of tectonic 
forcing F^. The lines denoted 'max' and 'min' indicate the boundaries of the limit cycles. Unstable 
fixed points or limit cycles are denoted by dashed lines. Calculations and figures were made using 
the pseudo- arc-length continuation software package AUTO [50]. (Online version in colour.) 



average concentration in CO2, consistently with an earlier proposal [49]. This 
tectonically driven decline is modelled here as a slow decrease in the forcing term 
F^{t) throughout the Pleistocene. 

Consider the bifurcation diagrams of the SM90 and SM91 models with respect 
to assuming no astronomical forcing (Fj = 0) (figure 3). The systems are then 
said to be free or autonomous. Depending on F M , both models show regimes with 
a stable fixed point, and regimes for which the fixed point is unstable, so that the 
system orbits along a limit cycle. 

These considerations led Saltzman to interpret the Middle Pleistocene 
Transition as a bifurcation between a 'quasi-linear' response regime to the 
astronomical forcing (in the fixed-point regime) to a regime of nonlinear 
synchronization (resonance) on the astronomical forcing. He concluded that ice 
ages would occur today even in the absence of astronomical forcing. The main 
effect of the astronomical forcing is to control the timing of glaciations. 

Saltzman's theory is seductive because it explains in a consistent framework 
several aspects of the Pleistocene climate history, including the change from 
linear to nonlinear regime [8], the presence of 100000 year periodicity in climate 
records [51], the lack of a 400 000 year spectral peak in the ice volume record 
(such a peak appears in the simple piece-wise linear model devised by Imbrie & 
Imbrie [52], owing to rectification of the precession signal), the synchronization 
of deglaciations on the astronomical forcing [53,54], and the occurrence of 
large climatic transitions even when eccentricity, which modulates the effect of 
precession on insolation, is at its lowest. 

The difficulty for accepting Saltzman's models as a definitive theory lies in 
the physical interpretation of the CO2 equation. This equation encapsulates all 
the interesting dynamics of the system and it is thus crucial to the theory. Some 
semi-empirical justification for the CO2 equation is given in Saltzman & Sutera 
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time (fraction of 76 ka) time (fraction of 1 10 ka) 

Figure 4. (a, b) Periodic and fixed-point solutions for SM90 (with = 1.2 ppmka ~ 1 ) and 
(c,d) SM91 (with = 0.5ppmka _1 ), near the sub-critical Hopf bifurcation points. The dynamics 
of SM90 slow down near the unstable fixed point, whereas the limit cycle of SM91 is much more 
decoupled from the position of the fixed point. (Online version in colour.) 



[44], but the form of this equation has undergone somewhat ad hoc adjustments 
in SM90. The form present in SM91 is again different, with important effects 
on the bifurcation structure, while the authors did not justify this latter change 
based on physical or biogeochemical considerations. 

To better appreciate the structural differences between the two models, 
let us return to the bifurcation diagrams. Consider SM90. As the forcing is 
decreased, the fixed point gives rise to a locally unstable limit cycle. The system 
must, therefore, find a stable limit cycle further away from the fixed point, but 
in this case not much further. This stable limit cycle is under the influence of 
the unstable fixed point and, in particular, the system slows down when it passes 
near it (figure 4). This scenario is depicted in figure 2 c. The limit cycle evolves 
as the tectonic forcing is further decreased, until it shrinks smoothly around a 
perpetually glaciated state. 

In SM91, the bifurcation induced by the decrease in tectonic forcing is much 
more explosive. The system lands on a stable limit cycle that turns out to be 
little affected by the position of the unstable point. The cycle dynamics do 
not show clear phases of acceleration and the system cannot be regarded as a 
relaxation system. The limit cycle disappears abruptly as the tectonic forcing is 
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Figure 5. Two histories of ice volume generated with the same models (a) SM90 and (b) SM91, 
using the astronomical forcing and a decline scenario for the tectonic forcing of CO2. The scenarios 
were chosen here to evidence the rise and decline of 100 000 year ice ages and are not the same 
as those used by Saltzman. Observe the explosive character of the appearance and decline of 
100 000 year cycles in SM91, with early and late excitations of the limit cycle by the astronomical 
forcing, and the smoother character of the evolution on SM90. Time is running from right 
to left. 



further decreased, through a phenomenon called a saddle bifurcation of cycles. 
The consequences of the difference between the bifurcation structures of SM90 
and SM91 may be further appreciated in the transient experiments shown 
in figure 5. 



(ii) Paillard's (1998) ice age model (P98) 

Paillard & Labeyrie [55] have been advocating the concept of relaxation for 
understanding palaeoclimate dynamics, both ice ages and the more abrupt events, 
since the publication of a seminal paper in 1994. We return to this article later 
on, and concentrate on another article published in 1998, in which Paillard [56] 
introduces a conceptual model of ice ages. Ice volume dynamics respond to an 
ordinary differential equation, 

dx = x R (y)-x _ ^ pgg ^ 
dt r R (y) T f 

In this equation, the ice volume x is linearly relaxed to xr with characteristic 
relaxation time tr. This relaxation process is further perturbed by the 
astronomical forcing F(t) with a characteristic time Tf. Such a system is said to 
be hybrid [57] because the relaxation equation involves a discrete state variable, 
here denoted by y. Its state may be 'deep glacial' (G), 'mild glacial' (g) or 
'interglacial' (z). The numerical values of and tr depend on this climate state. 
Climate states y follow a sequence i^ g —> G according to a set of conditions 
formulated on the level of glaciation x and insolation. Namely, the transition 
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g —> G is triggered when the forcing F{t) exceeds a certain threshold. Occurrence 
of G drives climate quickly into an interglacial state i because x&(G) and tr(G) 
are specified in the model to be low. 

Paillard is not very specific about the physical meaning of the discrete variable, 
but it accommodates the paradigm that the Atlantic Ocean circulation has 
gone through three different states during the latest glacial period: intermediate 
circulation, shut down of the circulation and modern, deep-sinking circulation. 
The system (P98) features the concept of slow-fast relaxation dynamics. However, 
this is not an oscillator because the shift from g to G is determined by the 
course of the external forcing. The Middle Pleistocene Transition is induced 
in (P98) in a fashion similar to that in Saltzman, and on the basis of similar 
physical assumptions (tectonically driven decline in CO2). The drift in climatic 
conditions induced by tectonics is accounted for by a term added to the 
astronomical forcing. In a later review, Paillard [58] further emphasizes empirical 
evidence for the relevance of the relaxation concept in the phenomenon of 
deglaciation. 

(hi) The Gildor-Tziperman model 

Gildor & Tziperman [59] take a moderate step towards higher model 
complexity by considering a slightly more explicit representation of atmosphere, 
ocean, sea ice and land ice dynamics. Namely, the ocean is divided into eight 
boxes, and the atmosphere into four boxes. The sea ice fraction responds to 
standard energy balance equations. More crucially, land ice growth is influenced 
by a somewhat controversial feedback between sea ice and precipitation. The 
feedback is controversial because it is assumed that cold climate results in a 
reduction in ice volume: sea ice growth causes a reduction in precipitation in ice- 
covered areas and, by this mechanism, almost suppresses accumulation of snow on 
ice sheets. The latter then no longer compensates for ice ablation and ice volume 
shrinks. 

A free oscillation arises from the fact that the ice volume thresholds for 
switching sea ice cover 'on' and 4 off differ. In other words, sea ice displays a 
hysteresis response to variations in ice volume. This is exactly the principle of 
the slow-fast relaxation oscillator depicted in figure 2 a: the curve of equilibrium of 
sea ice with respect to ice volume is the slow manifold, and ice volume integrates 
the state of sea ice in time. In turn, this oscillation can be synchronized on the 
astronomical forcing. 

The Gildor-Tziperman model is coupled to a biogeochemical cycle in a 
companion paper [60], but the essential dynamics of the glacial oscillation are 
unchanged. Tziperman et al [61] further comment on the model and its property 
of synchronization on the astronomical forcing, and find that its behaviour is 
essentially reducible to a hybrid dynamical system. 

(iv) The Paillard-Parrenin model 

Paillard & Parrenin [62] proposed yet another relaxation model in 2004 (PP04). 
The prognostic variables are ice volume /, the area of the Antarctic continental 
ice sheet A and the atmospheric concentration in CO2 (fi) (a . . . j are parameters): 
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dl 
dt 
dA 
dt 
d/x 
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l(dF(t) -e/ + /#(-£) + £- M ) 



1 



(-a/i- bF{t) + c-I), 



(PP04) 



and 



hi -iA+ j. 



As in the other ice age models, ice volume is a slow variable driven by 
the astronomical forcing. It is here coupled to a variable with a similar 
time scale (r^~ri) and a faster one (r /x = ri/3). The term H(—D), where H 
is the Heaviside function, represents the ventilation of the Southern Ocean. 
CO2 is released into the atmosphere when the Southern Ocean is ventilated 
(D<0), which drives deglaciation. Ice then grows slowly, until a Southern 
Ocean ventilation flush sends the system back to interglacial conditions. Ocean 
ventilation is thus the fast process in this model and it is the only nonlinear 
process accounted for. However, contrary to the Gildor-Tziperman model, it does 
not present a hysteresis behaviour. Consequently, the glacial cycles featured by 
this model cannot be interpreted in terms of shifts between the branches of a 
slow manifold. 

To better understand the dynamics of glacial cycles in this model, we consider 
the bifurcation diagram along typical solutions in the phase space for the free 
(i.e. unforced) system (figure 6). The parameter g is taken in this example as the 
control parameter, in order to preserve Saltzman's idea that ice age cycles appear 
as the consequence of a slow perturbation of the carbon cycle. As in SM90, PP04 
exhibits a limit cycle arising from a subcritical Hopf bifurcation. The dynamics 
along the limit cycle close to the bifurcation point are strongly influenced by 
the presence of the unstable focus. This is the configuration shown in figure 2c. 
Depending on g, the focus is either on the low-ice- volume side of the limit cycle 
(i.e. the system spends most of its time with high CO2) or on the high-volume side 
of the limit cycle (i.e. the system spends most of its time in low CO2). Parrenin 
and Paillard estimate that we are currently in the second configuration. 

(v) A minimal model of ice ages 

It has been claimed [61,63] that any model that has some form of 100000 year 
internal periodicity could be used to reproduce the course of ice volume over the 
last 800 000 years. Taking the argument at face value, Crucifix [64] used one of the 
simplest possible slow-fast oscillators — the van der Pol (VDP) oscillator — with 
minimal modifications to account for the astronomical forcing and the asymmetry 
between the phase of ice build-up and melt during the Late Pleistocene (a, /?, y 
and r are parameters; F(t) is the astronomical forcing), 



dx (-y + P + yF(t)) 



and 



dt r 

dy -a(y 3 /3 - y - x) 



(VDP) 
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Figure 6. Bifurcation diagrams and phase-space trajectories of the free Paillard-Parrenin model 
(a, b) PP04 and (c,d) the van der Pol (VDP) model. Both systems display explosive bifurcation 
scenarios, but the details are different. PP04 exhibits a subcritical Hopf bifurcation, whereas the 
limit cycle of VDP is explicitly framed by a slow manifold. Phase-space trajectories are drawn 
near the bifurcation points, that is: g = 0A in PP04 and (3 = 0.9 in VDP. The Heaviside function 
in PP04 is approximated as H(x) = a tan(500rc)/7r + 0.5 for analysis by AUTO. (Online version in 
colour.) 



The system dynamics are determined by the structure of the slow manifold 
x = y 3 /3 — y. The parameter (3 controls the position of the fixed point on the 
slow manifold and, consequently, the ratio of times spent by the system in the two 
branches ('glacial' and 'interglacial') of the slow manifold. The ice age curve can 
be captured with some tuning (figure 7), although it is fair to add that a small 
change in parameters may shift the timing of one or several ice age cycles. This 
minimal model was used to challenge intuitive arguments about the predictability 
of ice ages [64] . 

(b) Models for millennial climate variability 

(i) Dansgaard-Oeschger events as relaxation oscillations 

Welander [65] introduced the concept of relaxation oscillations in the context of 
ocean dynamics. He described a heat-salt oscillator involving exchanges of heat 
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Figure 7. Astronomical forcing, x and y trajectories obtained using system (van der Pol, VDP) with 
a = 30, /3 = 0.75, 7 = 0. 4 and r = 36 ka (1 ka = 1000 years). Black dots are an authoritative natural 
archive thought to mainly represent fluctuations in ice volume and deep-ocean temperature, and 
compiled in Lisiecki & Raymo [9]. 



and salt within a single oceanic column, coupled to a phenomenon of surface 
temperature relaxation. The destabilization process needed for the relaxation 
oscillation to appear is here related to diffusion between the deep ocean and the 
mixed layer. The system dynamics are further controlled by the mean freshwater 
flux at the top of the ocean column. It determines the transitions from a regime 
characterized by perpetual convection in the oceanic column, to a regime with 
intermittent convection (oscillation) and finally to a regime with no convection 
[66]. The bifurcations between the different regimes bear the character of global 
bifurcations, with the oscillation period approaching infinity near the bifurcation 
points (in particular, the second bifurcation bears the character of a homoclinic 
bifurcation). The heat-salt oscillator belongs thus to the class shown in figure 2b. 

Welander [67] and Winton & Sarachik [68] later introduce the concept of 
another kind of relaxation oscillator in the ocean. It involves the meridional 
structure of the ocean thermohaline circulation, and the key nonlinear process is 
the meridional advection of heat and salt. The oscillations featured by this model 
are termed 'deep-decoupling' oscillations [68]. Given that the slow process now 
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relates to heat accumulating in the global ocean, the characteristic time of deep- 
decoupling oscillations is of the order of 1000 years. The net flux of freshwater 
delivered to the North Atlantic acts as a bifurcation parameter controlling the 
transition between non-oscillating and oscillating regimes in the Winton-Sarachik 
model [69]. 

Millennial oscillations have since been observed across a hierarchy of ocean 
models, including three-dimensional ocean models with prescribed freshwater 
flux and restoring conditions to surface temperature [70] , and three-dimensional 
models coupled to a simple atmosphere [71,72]. Sakai & Peltier [73] proposed that 
millennial deep-decoupling oscillations could explain Dansgaard-Oeschger events. 
Colin de Verdiere et al. [33,74,75] complement this early proposal with a fairly 
complete theory based on ocean circulation model experiments. The oscillations 
described by Colin de Verdiere et al. involve the processes of turbulent vertical 
mixing (neglected in Winton & Sarachik [68]), advection and convection, which 
unify the salt oscillator with the deep-decoupling oscillation model. Incidentally, 
Colin de Verdiere et al. [74] dismiss the nonlinearity of the equation of state as 
the cause of the oscillations. 

There is, across the model hierarchy, consistency about the fact that the 
transition between the oscillating circulation regime and the so-called diffusive, 
haline regime (without deep convection) is associated with a homoclinic 
bifurcation [33,76]. The nature of the bifurcation between the convective regime 
and the oscillation is more model-dependent. Timmermann et al. [76], based on 
experiments with the eight-ocean box Gildor-Tziperman model, And a Hopf 
bifurcation; salt-conserving experiments with a two-dimensional ocean model 
show a transition towards a finite-period cycle, but of increasing period as the 
bifurcation is approached; experiments with a more idealized model, formulated 
as a two-equation dynamical system, reveal the signature of an infinite-period 
bifurcation [33]. 3 The latter implies that Dansgaard-Oeschger events, at the 
time when they appear soon after the glacial inception process, should be very 
long but of a similar amplitude to the Dansgaard-Oeschger events coming later 
in the glacial cycle. This feature is consistent with the Greenland ice core 
record (figure 1). More specifically, the first Dansgaard-Oeschger cycles that 
appeared at the beginning of the glacial era were characterized by a long 'plateau' 
phase (also called interstadial) during which the thermohaline circulation was 
certainly very active [15]. In the Colin de Verdiere et al. [74] theory, the 
plateau phase is the phase of the trajectory influenced by the 'ghost of the 
saddle point'. 

(ii) Dansgaard-Oeschger cycles as the manifestation of an excitable system 

Given the explosive nature of the bifurcations involved in ocean dynamics 
it is no surprise to find excitability properties in ocean models. Weaver & 
Hughes [70] discuss this effect in salt-conserving experiments with an 

3 This particular case is not illustrated in figure 2. There is no saddle point along or near the orbit, 
but there is a combination of parameters for which a fixed point appears on the limit cycle. Beyond 
this particular parameter value, this fixed point splits into a saddle and a node. This particular 
parameter value corresponds to an 'infinite-period' bifurcation. In practice, as long as the limit 
cycle exists, the trajectory slows down near the point where the saddle-node will appear. Some 
authors then refer to the influence of the 'ghost' of the saddle point [38]. 
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idealized-geometry, ocean model. The ocean-atmosphere model of intermediate 
complexity, CLIMBER (CLIMate BiosphERe model), was shown to exhibit 
excitability properties when boundary conditions are set to be typical of the latest 
glacial era [77]. The ocean circulation has then one stable state, with a moderate 
Atlantic overturning and a 'quasi-stable state' with more intense overturning. The 
conceptual sketch of the excitation cycles shown by Ganopolski & Rahmstorf [78, 
fig. 1] can be interpreted in terms of slow-fast dynamics, in which the different 
states of the ocean circulation constitute the different branches of a slow manifold. 
The intense overturning state, which is the 'plateau' phase of the Dansgaard- 
Oeschger event, may thus be viewed as the repelling branch of the slow manifold 
in the excitable regime (figure 2d). The excitable Dansgaard-Oeschger hypothesis 
was used as a possible basis to explain how a weak forcing, exogenous to the 
system, could explain the observed 1500 year periodicity of Dansgaard-Oeschger 
cycles (on this periodicity refer to earlier studies [79,80], but see the other view 
in the study of Ditlevsen & Ditlevsen [81]). Two such theories were developed on 
the basis of experiments with CLIMBER. One suggests that Dansgaard-Oeschger 
events are excited by stochastic fluctuations, modulated by a weak, hypothetical 
solar periodic forcing [82] (more on the effects of stochastic fluctuations in §4). 
The alternative theory suggests that the excitation is induced by the interference 
between two solar forcings with periods close to 1470/7(= 210) and 1470/17(~ 87) 
years [83], possibly combined with noise [84]. 

(hi) Heinrich cycles as a relaxation oscillation 

MacAyeal [85] proposed an ice binge-purge theory to explain Heinrich events. 
The theory rests on experiments with a one-spatial direction model of ice flow 
dynamics. Suppose, as a starting point, that ice volume grows in response to 
net accumulation of snow. The growth continues until the accumulated effect of 
geothermal heat flux causes basal sliding. A volume of ice is then released into 
the ocean (this is the 'purge'), causing the release of icebergs characteristic of 
Heinrich events. Ice volume thus decreases, until ice accumulation wins over so 
that ice volume can grow again. The ice binge-purge model is thus a relaxation 
oscillator combining a slow integrating process (ice mass accumulation) with a 
fast lateral discharge process. 

(iv) Coupling between Heinrich and Dansgaard-Oeschger events 

To what extent may Heinrich events interfere with Dansgaard-Oeschger 
dynamics? Paillard [55,86] investigated this question by coupling the MacAyeal 
ice model — but reduced to ordinary differential equations by Galerkin 
truncation — with a three-ocean box model. The coupling simply assumes that 
ice released into the ocean causes a net freshening of the surface of the North 
Atlantic that alters the deep-ocean circulation. Paillard realized that this coupling 
could lead to fairly non-intuitive and complex effects, such as the succession of 
Dansgaard-Oeschger events of decreasing amplitude between Heinrich events. 
This succession is known in the literature on palaeoclimate records as Bond cycles 
[18]. Paillard also found that the oscillations are aperiodic in this model under 
certain parameter configurations. 
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The issue is further explored in Schulz et al. [69], based on the Winton- 
Sarachik ocean model, and in Timmermann et al. [76], based on the slightly more 
sophisticated Gildor-Tziperman [59] ocean model. The objective was to study 
the response of deep-decoupling ocean oscillations to prescribed Heinrich cycles. 
Schulz et al. [69] noted that deep-decoupling oscillations could be synchronized 
on the Heinrich cycles. Timmermann et al. [76] then proposed, on the basis 
of numerical experiments with a fairly idealized model, that Heinrich events 
excite Dansgaard-Oeschger cycles because the variation in ice volume caused 
by a Heinrich event modifies slowly the amount of net freshwater released in the 
ocean. In turn, they suggested, Dansgaard-Oeschger may have a control on ice 
volume growth. This yields a two-way coupling between Dansgaard-Oeschger and 
Heinrich events. 

Experiments with more comprehensive models of the ocean-ice sheet- 
atmosphere system [87,88] generally support the idea that the different water 
and heat fluxes involved in the different phases of ice build-up and iceberg 
release are quantitatively sufficient to support a coupling between ice sheets 
and ocean circulation during the latest glacial era. However, it was also 
noted that 'three-dimensional thermomechanical ice-sheet models are unable to 
satisfactorily reproduce the binge-purge mechanism without an ad hoc basal 
parameterisation' [89]. 

To address this difficulty, a theory in which Dansgaard-Oeschger events trigger 
Heinrich events was recently proposed [89]. The ice shelf plays a key role, in 
blocking the ice stream flow from the ice sheet to the oceans. Heinrich events occur 
when this ice shelf is broken; for example, under the influence of ocean sub-surface 
warming associated with a Dansgaard-Oeschger event. The resulting model is 
a system displaying a slow ice build-up — the Heinrich release cycle excited by 
fluctuations in ocean sub-surface temperature. 

(v) Holocene oscillations and relationship with Dansgaard-Oeschger events 

The much smaller ocean oscillations that characterized the Holocene period 
may also be a relaxation phenomenon. Schulz et al. [31] observed oscillations in 
the atmosphere-ocean model of intermediate complexity ECBILT-CLIO. These 
oscillations are related to the convective activity in the Labrador Sea. 

Schulz et al. [90] considered the existence of such an oscillator in an earlier 
reference and speculated on the possible interactions between the centennial 
oscillations, millennial oscillations and Heinrich cycles. They considered a model 
in which each of these three kinds of oscillations is modelled as a Morris-Lecar 
relaxation oscillator. Their working hypothesis is that glacial conditions induce 
a coupling between these oscillators. They then observed that a very stable 
1500 year oscillation appears, which they interpreted as a model equivalent of 
Dansgaard-Oeschger events. 

4. Stochastic effects 

The myriad of chaotic motions that characterize the dynamics of the ocean 
and the atmosphere may be taken into account in the form of parametrizations 
involving stochastic time processes. The method was introduced in climatology 
in the 1970s [91] and the theoretical justifications, which allow one to model 
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chaotic motions as a (linear) stochastic process, are reviewed by Penland [92,93]. 
In a statistical inferential framework, the stochastic parametrizations may also 
be viewed as a way to account for the distance necessarily existing between 
the concepts and dynamics featured by the model, and the complex system 
being observed. 

The effects of stochastic processes on relaxation oscillators and excitable 
systems are generally well documented in the literature because this is a topic 
of general interest [94]. Here, we review some of them in the specific context of 
palaeoclimate dynamics. 

(a) Stochastic effects on ice age dynamics 

(i) Phase dispersion 

One of the basic effects of noise on oscillators is the phenomenon of phase 
dispersion: a weak stochastic forcing on an oscillator causes a fading out of 
the memory of the exact initial conditions, even though the gross structure of 
the oscillation visualized in the phase space is conserved. The phenomenon is 
well known and it is an immediate consequence of the neutral stability of the 
phase of a free oscillator with respect to fluctuations. It was previously suggested 
that this phenomenon of phase dispersion may concern ice ages [95], but it is 
more commonly believed that ice ages are phase locked on the astronomical 
forcing. This phase locking should act against dispersion and permit a very long 
predictability horizon of ice ages. However, a phenomenon of phase dispersion may 
happen in oscillators that are locked on a periodic forcing. A stochastic fluctuation 
may momentarily cause a burst of desynchronization, called phase slip, during 
which the system is unhooked from its corresponding deterministic trajectory and 
attracted to another trajectory, which leads or lags the original one by one forcing 
period [41, §3.1.3]. The difference between phase diffusion in a free system and in 
a periodic forcing-driven oscillator is that the diffusion effect has, in the latter, a 
quantum nature. More formally, it is said that the stochastic forcing disperses the 
system states around the different attractors that are compatible with the forcing. 
In a work in preparation, we suggest that the astronomically forced climate system 
may satisfy the conditions for a similar phenomenon of phase dispersion to occur 
(B. De Saedeleer, M. Crucifix & S. Wieczorek 2011, unpublished data). Given that 
the astronomical forcing is aperiodic the description of the phenomenon requires 
a suitable theoretical framework, which relies on the notion of a 'local pullback 
attractor'. The equivalent of a phase slip is, in the aperiodic forcing context, 
a stochastic shift from one of the deterministic pullback attractors to another 
one. The phenomenon is illustrated based on experiments with the VDP model 
in figure 8. 

(ii) Reduction of period 

Additive fluctuations generally reduce the period of relaxation oscillators. In 
an oscillator presenting a homoclinic orbit such as the Duffing oscillator, additive 
fluctuations reduce the time spent near the unstable focus [96]. This implies that, 
even at the corresponding bifurcation point in the deterministic system, the return 
time of oscillations in the stochastically perturbed system remains finite. In a 
slow-fast oscillator such as the VDP oscillator, additive fluctuations generally 
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Figure 8. The phenomenon of trajectory shift illustrated with the model (van der Pol, VDP) with 
astronomical forcing, with parameters r = 36, = 0.75, y = 0.4, a = 30 ka. The black line represents 
the deterministic system. The generated history reproduces reasonably well the fluctuations of ice 
volume of the last 800 000 years. The red (grey) line represents the same system but with an 
additive Wiener process added to the fast variable y, with variance b = 0.2 (ka)" 1 / 2 . One possible 
realization of the stochastic system is shown. Observe the solution shift around 450 000 years ago. 
Depending on the realization, the shift may occur at different places. (Online version in colour.) 

result in early escapes of the branch of the slow manifold on which the system 
lies (figure 9). The period of the oscillator is thus affected by a correction that 
increases approximately linearly with the noise variance in the slow-fast VDP 
oscillator [97]. 4 This property was used at least once in Pleistocene theory, in the 
silicate weathering hypothesis advanced by Toggweiler [99] . Additive fluctuations 
reduce the limit cycle period from 800000 to about 100 000 years. The reasons 
for the period reduction being so dramatic are left for another study. 

4 This result is established assuming fluctuations added to the slow variable. A much more general 
theory, suitable for different slow manifolds and additive fluctuations to slow and fast variables, is 
now available [98]. 
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Figure 9. Two possible effects of additive noisy fluctuations in a slow-fast system, (a) If the 
slow-fast system is in the excitable regime, then noisy fluctuations can cause excitations. 
Excitations will be sporadic if the amplitude of the oscillations is weak (as shown here), or regular, 
as in a limit cycle, in the coherent resonance regime, (b) If the slow-fast system is oscillating, 
then additive fluctuations may cause early or delayed escapes from the slow branches, compared 
with the deterministic system. The net effect is a shortening of the cycle duration. (Online version 
in colour.) 



(b) Stochastic effects on Dansgaard-Oeschger dynamics 

(i) Stochastic excitation and resonance 

Noise may naturally act as an excitation agent in an excitable system (figure 9). 
The topic is extensively reviewed by Lindner et al. [94]. Excitation loops are 
sporadic if the noise amplitude is weak, in which case the recurrence time of the 
events is set by the noise amplitude, whereas the amplitude of the events is set by 
the structure of the deterministic vector held. The frequency of excitation loops 
increases with the noise amplitude, until the system behaviour is qualitatively 
similar to a limit cycle regime. This is the coherent resonance regime. For yet 
higher noise amplitudes, the limit cycle structure is destroyed. 
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The concept of stochastic excitation has been considered several times in 
Dansgaard-Oeschger theories. The idea is introduced based on experiments with 
an ocean general circulation model with idealized geometry and forcing [70] . The 
effect of stochastic fluctuations is not only to excite oscillations in a system 
normally at rest, but also to reduce the period of these oscillations when the 
system is in oscillatory regime. 

A phenomenon of non-autonomous stochastic resonance may occur if the noise 
is superimposed on a weak external drive. For this to happen, the autonomous 
system needs to be stable but excitable. The external forcing must be too weak to 
cause excitation by itself. The role of the noise is to provide the additional power 
to induce excitation. The timing of the excitation is then related to the phase of 
the external forcing. The mechanism was proposed several times [76,82] to explain 
the 1500 year recurrence time of Dansgaard-Oeschger events. The idea remains 
questioned, either on the grounds that the 1500 year recurrence time observed in 
palaeoclimate records is coincidental [81] or on the grounds that the 1500 year 
external forcing is unidentified [33,83]. A more subtle case of stochastic resonance 
involves the combination of noise with two solar cycles of 210 and 87 years, which 
yields the concept of 'ghost resonance' [100] for which some support, albeit not 
conclusive, is found in the observations [80]. 

(ii) Decreased sensitivity to noise in resonant oscillators 

Coupled oscillators may exhibit, collectively, a resonance period that is more 
robust to external fluctuations than the uncoupled oscillators. Schulz et al. [90] 
used this property to explain the stability of the Dansgaard-Oeschger recurrence 
period of 1500 years in the presence of random fluctuations, without having to 
invoke an external forcing. 

(hi) Pseudo-oscillations in two-well systems 

Finally, a behaviour reminiscent of oscillations may occur in a system that 
is neither oscillating nor excitable, but which presents several stable states. 
Noise then simply induces jumps between these different states. The simplest 
mathematical model is the Langevin equation and this is on this basis that Schulz 
et al. [31] interpret the Holocene oscillations observed in the ECBILT-CLIO 
climate model. 

5. Concluding discussion: can dynamical systems be used for inference? 

The review has shown that relaxation oscillations are a popular and powerful 
model to explain oscillations observed in the Pleistocene record. The concept of 
relaxation implies some form of slow-fast separation, in the sense that at least 
one component of the system spends most of its time in 'quasi-equilibrium' states 
(this may be a 'slow manifold branch' or a region influenced by a saddle point, 
depending on the system structure), with acceleration phases. 

Some of these models were constructed following a fairly careful procedure 
of truncation of a system of partial differential equations, which describes 
some of the fluid dynamics of the climate system. Others were proposed on a 
more conceptual basis, the idea being precisely to test a hypothesis based on 
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palaeoclimate observations. The latter approach is sometimes criticized, on the 
grounds that box models, for example, cannot reasonably be taken as an adequate 
representation of the complex dynamics of the oceans [101]. 

This leads us to the last question of this review: can dynamical systems be 
used for inference on palaeoclimates? Inference implies that something is being 
learned by confronting a model with observations. This inference process may 
take the form of a calibration procedure (update our knowledge on parameters 
on the basis of observations) or a model selection procedure (which model, among 
different alternatives, explains the observations best). 

The position taken here is that there is not such a thing as an 'attractor' 
of the climate system that is to be 'discovered'. The hope is that some of its 
modes of behaviour are sufficiently decoupled from the rest of the variability 
to justify the fact that simple dynamical systems may capture the fundamental 
dynamical properties of these modes, and we want to learn about these modes 
from palaeoclimate observations. 

The programme is challenging. Indeed, it was emphasized that different 
physical assumptions may lead to dynamical systems with dynamical properties 
that are similar enough to produce a convincing visual fit on palaeoclimate data 
[61]. The message is largely echoed in this review. The modeller's challenge is, 
therefore, to operate a model selection on more stringent criteria than just fitting 
some standard time series. For example, palaeoclimate observations may yield 
constraints on the bifurcation structure of the system. The Middle Pleistocene 
Transition is an attractive test case in this respect. 

In a statistical inference process, the observations should be a plausible outcome 
or realization of the model. This makes sense only if the model has a stochastic 
component, which describes its uncertainties, limitations, and the noise that 
emerges from the chaotic motions of the atmosphere and oceans. 

Stochastic dynamical systems are beginning to be used for inference on 
palaeoclimate time series. In a method called 'potential analysis', the climate 
system is modelled as a Langevin equation, that is, the combination of a down- 
gradient drift with a Wiener additive process and inference is made on the number 
of wells of the potential function [102]. The method was applied on Pleistocene 
climate records, yielding the conclusion that the number of wells increased from 
two to three over the course of the Pleistocene [103]. 

However, our position so far has been to favour a Bayesian methodology, 
because it allows one to encode physical constraints in the form of prior 
distributions on model parameters. The Bayesian formalism is also naturally 
designed for model calibration, selection and probabilistic predictions. 

The fact is that Bayesian methods for selection and calibration of dynamical 
systems on noisy observations are only emerging. In a recent attempt, we 
considered a particle filter for parameter and state estimation [104]. To be honest, 
there is ample room for progress. Whether the process of inference with simple 
dynamical systems on palaeoclimate data will lead new insight in this context 
still needs to be demonstrated. 

This review benefited from stimulating discussions during the Isaac Newton Programme 
Mathematical and Statistical Approaches to Climate Modelling and Prediction held in Cambridge 
during the summer-autumn 2010. Guillaume Lenoir, Bernard De Saedeleer and Jonathan Rougier 
provided helpful comments on a first version of the article. Thanks are also due to Didier 
Paillard, Olivier Arzel, Alain Colin de Verdiere, Andre Paul and Sebastian Wieczorek for e-mail 
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